Observation of non-sinusoidal current-phase relation in graphene Josephson junctions 
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We present direct measurements of the current-phase relation for lateral Josephson junctions 
with a graphene barrier, obtained by a phase-sensitive SQUID interferometry technique. We find 
that the current-phase relation is forward skewed with respect to the commonly observed sinusoidal 
behavior for short junctions in the quasi-ballistic transport regime, consistent with predictions for 
the behavior of Dirac fermions in a Josephson junction. The skewness increases with critical current 
and decreases sharply with increasing temperature. 



The interplay of superconductivity and the unique 
electronic structure of graphene leads to unusual co- 
herence effects such as gate-tunable supercurrents [1], 
a finite supercurrent at the Dirac point [2], and spec- 
ular Andreev reflection [3]. Much recent work has fo- 
cused on superconductor- graphene-superconductor (S-g- 
S) Josephson junctions in which theoretical [4, 5] and 
experimental studies [6-8] have examined the effects of 
parameters such as the junction geometry and number of 
layers on the critical current. However, unique informa- 
tion about the processes influencing the supercurrent can 
be obtained by measuring not just the magnitude of the 
supercurrent but also its dependence on the phase dif- 
ference across the junction, characterized by the Joseph- 
son current-phase relation (CPR). The simplest models 
of Josephson tunneling predict a sinusoidal variation of 
the current with phase; deviations such as skewness have 
been predicted [9], but rarely observed. It is possible to 
extract some information about the CPR of a junction by 
measuring critical current diffraction patterns, Shapiro 
steps, or switching currents in a SQUID (superconduct- 
ing quantum interference device) configuration [10], but 
the most definitive approach is to measure the CPR di- 
rectly using a phase-sensitive interferometry technique. 

In this Letter, we present experimental measurements 
of the CPR in Josephson junctions having a single-layer 
graphene barrier. The junction is incorporated into a 
superconducting loop coupled to a dc SQUID which al- 
lows the junction phase to be extracted directly [11] (see 
Fig. 1(a)). We observe significant deviations from the typ- 
ical sinusoidal behavior for short S-g-S junctions, where 
the junction length {L) is less than the superconduct- 
ing coherence length (^) in the junction. The deviations 
consist of a forward (positive) skewness in the CPR that 
varies as a function of critical current and is detectable 
at temperatures below 300mK. This behavior is different 
from typical metallic Josephson junctions [9, 12], and ap- 
pears to be well-described by self-consistent tight-binding 
Bogoliubov de-Gennes (TB-BdG) calculations which con- 
sider the Dirac spectrum of graphene [13]. The CPR 
curves of junctions with L > ^ do not exhibit significant 
skewness. 

Measurements were performed on four junctions on 
four separate graphene flakes. The sample dimensions 
and characteristics are given in Table 1. All samples were 
prepared by mechanical exfoliation [14] of graphite flakes 




(a) Interferometer Setup (b) S-g-S 



FIG. 1. (a) Circuit diagram of the current-phase interferome- 
try experiment. The graphene is depicted by the honeycomb 
lattice and the junction is shown on top. Note that the in- 
ductance is coupled to the SQUID through a flux transformer 
(not shown), (b) SEM image of the S-g-S junction. Labels 
SI, S2, CI, C2 correspond to identical labels in figure (a). 



Sample 


L(nm) 


Ty(/im) 




/co (nA) 


Jco(nA//im) 


A 


70 


0.3 


0.35 


71 


237 


B 


100 


0.5 


0.5 


107 


213 


C 


350 


3 


1.75 


39 


13 


D 


350 


10 


1.75 


160 


16 



TABLE I. Table summarizing the characteristics of each junc- 
tion. I CO and J CO are the critical current and critical current 
density {Ico/W) of the junctions at the Dirac point at low 
temperature (T < 50mK). 



onto highly p-doped Si substrates covered by 300nm of 
Si02, where the doped substrate acts as a global back- 
gate. Samples were annealed at a temperature of 400C 
in 1900/1700 seem H2/Ar, and subsequently character- 
ized by optical imaging, atomic force microscopy, and Ra- 
man spectroscopy. Ti(4nm)/Al(60nm) contacts were fab- 
ricated via electron beam lithography and electron beam 
evaporation. Samples A and B were fabricated on large 
flakes of graphene, with the area of the graphene much 
larger than the junction size (see Fig. 1(b) and Sup- 
plementary Information (SI)). Junctions C and D were 
fabricated on narrow strips of graphene, with the edges 
of the graphene naturally defining the width of the junc- 
tions. Measurements were performed in a dilution refrig- 
erator having a base temperature of lOmK. Transport 
properties of the graphene flakes were estimated from the 
normal state conductivity, which gave mean free paths 
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FIG. 2. The current phase relation (CPR) at (from top to 
bottom) lOmK, 150mK, and 350mK in the short junction 
hmit (L < ^, (a)) and the long junction limit (L > <J, (b)). 
The measured CPR curves are the solid lines. The dotted lines 
are theoretical curves taken from [15]. The scatter points are 
theoretically calculated in [13]. All theoretical data is scaled 
to fit the measured curves. 



(/) of ~ 60nm and mobilities (/i) of ~ 6000 cm^V~^s~^ 
(see SI). Samples A and B, where the junction length L 
is of order are considered to be in the quasi-ballistic 
regime, whereas samples C and D, where L ^ are con- 
sidered to be in the diffusive regime. The superconduct- 
ing coherence length in the junction is estimated to be 
i - .JhBjK - 200 nm, where D is the diffusion length, 
Vj \^ the Fermi velocity, and A is the superconducting 
gap of the electrodes. We can then characterize samples 
A and B, where L < ^, as in the short junction limit, and 
samples C and D, where L > ^, as in the long junction 
limit. 

The circuit used to extract the CPR is shown in Fig. 
1(a). The Josephson junction is connected in parallel 
with a fabricated thin film superconducting loop of in- 
ductance L. A current / injected into the circuit divides 
so that the phases across the junction and the loop are 



equal (i.e. ^ = 27r(^/^o)) where ^ is the flux in the loop 
and is a flux quantum. The component of current 
in the loop inductor is measured by coupling the flux in 
the loop to a commercial dc SQUID via a filamentary 
superconducting flux transformer. From the circuit in 
Fig. 1(a), it can be seen that the current Ij and phase (\) 
across the junction are given by: 



• = 27r 



y SQUID 



VSQUID 



(1) 



(2) 



where II is the current through the inductor, and 
VsQUiD is the measured SQUID voltage, with the flux 
transfer function {dVsquiD /d^)- To extract the CPR, 
it is necessary that the Vsquid vs / characteristics re- 
main non-hysteretic. For a sinusoidal CPR, this requires 
that Pl = 27tLIc/^o < 1; for a CPR with higher har- 
monics, the condition is more stringent, however due to 
noise-rounding this condition gives a reasonable estimate 
for the range of critical currents at which the CPR can 
be determined, as discussed in the SI. For our circuit, 
L ~2.3nH, which sets the limit of Ic < 135nA for the 
junction. 

CPR curves measured at the Dirac point for Sample A 
(L < ^) are shown in Fig. 2(a). Clear non-harmonic be- 
havior (forward tilt/skewness) is observed at T = lOmK 
and T = 150mK (similar behavior is observed in Sam- 
ple B). Also shown are curves and individual points cal- 
culated by the Dirac Bogoliubov de Gennes (DBdG) 
method and the tight binding Bogoliubov de Gennes 
(TB-BdG) method, respectively [2, 13] (these methods 
will be discussed further in the theory section). At all 
temperatures there is a good fit between the measured 
CPR curves and the CPR curves calculated by the TB- 
BdG method. The DBdG method provides a worse fit at 
low temperature, but converges to the measured result 
above 300mK. Measured curves from Sample C (L > ^), 
shown in Fig. 2(b), exhibit typical sinusoidal behavior 
at all temperatures (similar behavior is observed in Sam- 
ple D). Thus only the short junctions seem to exhibit 
non-sinusoidal, forward- skewed behavior. 

From the CPR measurements, the Ic can be extracted 
as a function of the back-gate voltage Vg- Fig. 3(a) 
displays the Vg dependence of Ic for sample A, along 
with the Vg dependence of the junction normal resis- 
tance i^AT at low temperature. The Dirac point resides 
at Vg = —35V (similar to sample B). We consistently ob- 
serve an asymmetry in Ic{Vg), common to all samples, 
that becomes more pronounced at higher temperatures, 
while Rn{Vg) appears to be symmetric near the Dirac 
point. Ic{Vg) increases sharply on the n-doped side of 
the Dirac point (Vg > —35V) , but increases slowly up to 
a constant value on the p-doped side {Vg < —35V). This 
behavior is not well-understood, but has been observed 
elsewhere [1]. In order to satisfy the constraint < 1 
and accurately measure the CPR, the Ic of the junction 
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should be below the dashed line shown in Fig. 3(a). This 
is achieved at lOmK and lOOmK for Vg < —20V, and 
is easily satisfied for higher temperatures. The Jc{T) 
dependence, where Jc = Ic/W^ is shown in Fig. 3(b) 
for all samples, along with Jc{T) curves calculated by 
the TB-BdG method in [13]. The significant change in 
Jc between junctions with L < ^ and junctions with 
L > ^ matches the theoretical prediction for short and 
long junctions. 

We parameterize the skewness (non-harmonic behav- 
ior) of the measured CPR curves by a variable S = 
{2(j)max/'7T) — 1, where (j)max IS the position of the maxi- 
mum of the CPR; S ranges from to 1 as the CPR evolves 
from a sine wave towards a forward saw-tooth wave. The 
Ic dependence of S is shown in Fig. 4(a) for the n-doped 
side of the Dirac cone for Sample A (Sample B shows 
similar behavior). Since Ic{Vg) in the p-doped regime is 
not well-understood, we omit data in this regime for sim- 
plicity (see SI for details). Above Ic = 135nA, Pl > 1 
and hysteretic switching behavior averaged by noise in 
the measurement circuit causes the CPR to appear neg- 
atively skewed, resulting in a sharp decrease in S vs Ic 
[see SI] . Note that this decrease is an artifact of the mea- 
surement itself, not a result of physical processes in the 
junction. However, below Ic = 135nA, < 1 allows 
for an accurate extraction of the CPR from the measure- 
ment. A positive skewness is observed in this regime that 
increases linearly from the Dirac point to Ic = 135nA 
for temperatures below 300mK. Above 300mK, a much 
weaker dependence of S on Ic is observed. The skewness 
also exhibits a strong temperature dependence (see Fig. 
4(b)). As T ^ Tc^ S approaches for all samples. S in- 
creases significantly below 300mK for samples A and B, 
while it remains close to zero for sample C. S for sample 
D (not shown) exhibits behavior similar to that of sam- 
ple C, but is negative below 250mK due to large critical 
currents {Ic > 135nA, > 1). 

The CPR of a ballistic S-g-S junction, at zero tempera- 
ture and in the short-junction regime, was first predicted 
theoretically by applying the standard Bogoliubov-de 
Gennes theory to the Dirac spectrum, also known as the 
DBdG formalism [2]. Ic is carried by Andreev bound 
states in the junction according to: 
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FIG. 3. (a) Critical current measurements (left axis) as a 
function of the back-gate voltage shown at lOmK up to 800mK 
for Sample A. The normal resistance at 1.2K (right axis), 
shown as the solid black line, is also given as a function 
of gate voltage. The dotted line at Ic =135nA (/^l = 1) 
marks the maximum critical current at which the CPR can 
be determined, (b) Measured critical current as a function of 
temperature for each sample. Solid lines are measurements, 
while dashed lines are theoretical fits taken from ref. 10. 
Jco=237nA//im is taken to be the Jc of sample A at lOmK. 
Jc is defined to be IcjW . 



ic=ff: (3) 

i=o ^l-TnSm\(l)/2) 

where A is the superconducting energy gap and are 
the transmission coefficients, which are functions of the 
Andreev bound state wavevectors (kn)- Large values of 
Tn lead to forward- skewed contributions to /c(^)- Near 
the Dirac point, propagating states (real kn) and evanes- 
cent states (imaginary kn) both contribute significantly 
to the conduction. There are only a few number of states 
with high T^, allowing states with lower Tn to reduce the 
skewness. Far from the Dirac point, many propagating 
states with high Tn dominate. A closed form solution to 
Eq.(3) at the Dirac point for widths W ^ L is given by: 



IM = ^^cos((/)/2)tanh-^(sin((/)/2)) (4) 

which has a forward skewness of 0.255 [2]. In the 
DBdG approach, when finite gate voltages are applied to 
the graphene, the number of propagating bound states 
increases and dominates the supercurrent. This causes 
the skewness in the CPR to first increase, then oscillate 
due to interference effects between the bound states, and 
then saturate at 5* ~ 0.42 with increasing carrier density. 
The DBdG approach assumes A is fixed at the leads of 
the superconductor. This boundary condition couples 
electron and hole transport in the junction, leading to 
Andreev reflection. 
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FIG. 4. (a) Measured skewness (S) of the CPR as a function 
of critical current (/c) for sample A. The horizontal dotted 
lines represent predictions for the skewness of the CPR at the 
Dirac point by the DBdG formalism and the TB-BdG for- 
malism. The vertical dotted line indicates the critical current 
at which /^l = 1. The solid and dashed lines are guides to 
the eye for points measured in the non-hysteretic {/3l < 1) 
and hysteretic regimes {13 l > 1), respectively, (b) Measured 
skewness as a function of temperature for all samples at the 
Dirac point. The black dashed line is a theoretical fit taken 
from [13] for the case of short junctions (L < ^) at the Dirac 
point. 



We observe skewness values up to 0.17, somewhat be- 
low the predictions of DBdG formalism. The DBdG ap- 
proach has been extended to arbitrary temperatures, but 
the results still predict skewness values higher than we 
observe [15]. This discrepancy may be due to the fact 
that the DBdG approach takes into account neither cur- 
rent depairing nor proximity effect in the superconduct- 
ing electrodes, both of which reduce A in the electrodes 
and have the effect of reducing the skewness. A different 
approach, previously referred to as the TB-BdG formal- 
ism, incorporates this effect by allowing A to vary in the 



electrodes as well as in the junction [13, 16, 17]. The re- 
sults from [13] indicate that current depairing affects the 
CPR predominantly in short junctions, while the proxim- 
ity effect affects the CPR in long junctions. We find that 
our short junction CPR data at low temperature fits rea- 
sonably well to CPR curves calculated by the TB-BdG 
method (see Fig. 2(a)) indicating that current depair- 
ing is likely a significant physical mechanism affecting 
the junctions. Calculations that include a temperature 
dependence, but not current depairing, are also shown 
[15]. However, these curves poorly match the data. It is 
also clear from Fig. 4(a) that typical measured skewness 
values at low temperatures are comparable to skewness 
values predicted by the TB-BdG approach (lower dashed 
line), while they are below those predicted by the DBdG 
approach (upper dashed line). Some of this difference 
may be attributed noise-rounding, which could decrease 
the measured skewness (see SI). 

Measurements of 5* vs T and a prediction by the TB- 
BdG formalism at the Dirac point are shown in Fig. 4(b) 
[13]. The temperature dependence for Sample A, mea- 
sured near the Dirac point, generally follows the theoreti- 
cal prediction, implying that the TB-BdG predictions are 
reasonable for our data. Sample B also qualitatively fits 
the prediction, though it exhibits slightly lower skewness 
values which we attribute to the longer L. The TB-BdG 
approach predicts negative skewness for short junctions 
at higher temperatures, particularly for large |Vg|. We 
observe only minor negative skewness values for samples 
A and B at higher temperatures. This is probably due 
to the restriction that Ic < 135nA, which limits \Vg\- 
The CPR for long junctions can also be calculated us- 
ing the TB-BdG approach, but our long junction CPR 
curves (samples C, D) do not match those calculated in 
[13] (see Fig. 2(b)) for low temperature (T < 300mK). 
However, it is uncertain whether the TB-BdG approach 
is applicable to junctions with diffusive {L I) trans- 
port. 

In conclusion, the CPR curves of graphene Josephson 
junctions with varying lengths have been measured via a 
phase- sensitive interferometry technique. Positive skew- 
ness values are reported for short junctions at low tem- 
peratures. The measured CPR curves of the short junc- 
tions are consistent with those calculated with the TB- 
BdG formalism. The results suggest that the CPR of a 
ballistic graphene junction is dominated by a few num- 
ber of propagating Andreev bound states. This behavior, 
which is characteristic of Dirac fermions, is different from 
typical metallic Josephson junctions. However, to pro- 
vide a complete picture of the dynamics of the junction, 
gap suppression due to current depairing in the electrodes 
should be accounted for. 
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FIG. SI. Optical image (80x) of sample B. The device in 
the image is covered with PMMA, with the leads defined by 
e-beam lithography. 



I. TRANSPORT PROPERTIES 

Transport properties such as / and ji are calculated 
from four-probe measurements of RiVc) on separate 
graphene flakes. We estimate / from / = cFh/2e^kF^ and 
11 from ji = aTr/ekp using kp = vT^^~^^^^o)^r£o7r7ed. 
Here d is the thickness of the Si02 substrate, a is the 
conductivity, Er is the dielectric constant, and kp is the 
Fermi wavelength. We use values of |Vg| > lOV from the 
Dirac point, where 1{Vg) and /i{Vg) are approximately 
constant and charge impurities affecting the conductivity 
near the Dirac point are screened [1, 2]. 



II. MEASUREMENT 

The circuit in Fig. 1(a) is suitable for a DC mea- 
surement, where the applied current / is increased lin- 
early while the output voltage of the SQUID (Vsquid) 
is recorded directly, yielding the dashed line in Fig. S2. 
Vsquid converts (via the flux transfer function V^, and 
inductance L) to the phase of the junction. To extract 
the CPR, the vertical and horizontal axes in fig. S2 
are inverted and a linear fit is subsequently subtracted, 
yielding the CPR. The DC measurement is highly sus- 
ceptible to flux noise and drift from the SQUID. We 
therefore choose to measure the differential of Vsquid^ 
or AVsQuiD^ via a lock-in technique. In this AC mea- 
surement, / is linearly increased while a small oscillation 
current (A/) from a lock-in amplifier sums with / to pro- 
duce a differential SQUID response AVsquid (fig- S2, 
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FIG. S2. Measurements obtained from the experimental 
setup in Fig. la. AVsquid (solid curve) is the differential 
SQUID voltage (AC coupled) measured by the lock-in am- 
plifier. Vsquid (dashed curve) is the SQUID voltage output 
(DC coupled). 



solid curve) that is detected by the lock-in amplifier in- 
put. We subsequently integrate AVsquid to produce 
VsQUiD' This approach filters the noise and drift men- 
tioned previously, allowing for a more precise CPR ex- 
traction. 



III. ASYMMETRY: P-DOPED VS. 

REGIMES 



AT-DOPED 



As noted in the letter, we only present skewness val- 
ues on the n-doped side of the Dirac cone. The com- 
plete set of skewness measurements (p-doped and n- 
doped regimes) is shown in figure S3 for samples A and 
B. Note that the skewness linearly increases on the n- 
doped side of the Dirac cone, while the skewness slightly 
decreases and then flattens out on the p-doped side. Cal- 
culations from [3] indicate that S{Vg) should be symmet- 
rical around the Dirac point, contrary to our observation. 
Asymmetries in p-type and n-type conductivities around 
the Dirac point in graphene have been explained by the 
difference in scattering cross-section between holes and 
electrons on charge impurities [4, 5], and also by p-p and 
p-n junctions formed by the leads at the graphene in- 
terface [6]. In ref. [6], Ti-Al electrodes are reported to 
produce a higher n-type conductivity than p-type con- 
ductivity, consistent with our observations. This suggests 
that the contacts are primarily causing the asymmetries 
observed in our measurements. 
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Backgate Voltage (V) 

FIG. S3. Skewness vs. backgate voltage for samples A and 
B measured at lOmK. The vertical dashed line indicates the 
location of the Dirac point. The solid lines are guides to the 
eye. 



IV. HYSTERESIS / NOISE ROUNDING 

In Figure 4(a), our results show a strong suppression of 
the forward skewness as the critical current is increased 
above Ic = 135nA, resulting in a backward-skewed OPR. 
We understand this effect as being due to noise-rounding 
in the CPR measurement. The phase interferometer 
technique depends on measuring the fraction of the ap- 
plied current flowing through the SQUID loop. This cre- 
ates the step pattern shown in Figure S2. For a sinusoidal 
CPR, these steps become sharp as /^^ approaches 1 and 
become hysteretic for /3l > 1. As a result, noise induced 
in the circuit can round out these characteristics, distort- 
ing the extracted CPR and even removing the hysteresis 
so that CPR curves can be extracted even in the regime 
for 13 L > 1. Noise effects can be substantial in the junc- 
tions because of the small critical currents. 

In Figure S4(a), we demonstrate this effect on a pure 
sin(^) CPR by modeling the SQUID current in the pres- 
ence of Gaussian noise in the applied current. The 
RMS current fluctuations {SI) of the noise are set to 
SI = 0.1(^0/^)7 a relatively large value chosen to exag- 
gerate the noise-rounding effect. The noise preferentially 
smears out the sharp portions of the steps, resulting in 
the backward-skewed CPR in Figure S4(b). This effect 
is even more pronounced for the intrinsically forward- 
skewed characteristics that describe the graphene junc- 
tions. For the predicted (and observed) CPR, the SQUID 
current curves are hysteretic for all values of so the 
noise affects the response in all regimes. An example of 
how this affects the CPR curve predicted by Equation 4 
is shown in Figure S4(c). In Figure S5, we use this noise- 



rounding model to show the effect of the level of noise 
on the CPR skewness as a function of {Pl Ic). As 
expected, the noise-rounding generates a backward skew- 
ness that suppresses and then dominates the CPR as /3l 
increases. Comparison with Figure 4(a) verifies that this 
mechanism explains this prominent feature of our data. 
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FIG. S4. (a) Modeled step response of Vsquid vs applied cur- 
rent for a sin(0) CPR. (b) Extracted CPR from (a), (c) Mod- 
eled step response of Vsquid vs applied current for Eq.(4). 
(d) Extracted CPR from (c). For (a)-(d), the solid line is the 
intrinsic response and the dashed line is the noise-rounded 
response. 
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FIG. S5. Skewness vs /Sl for a sinusoidal CPR and for the 
skewed CPR given by Eq.(4). 
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